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I.  INTRODUCTION 

The  forced  liquid  project  is  being  directed  by  the  Munitions  Branch  of  the 
U.S.  Army  Chemical  Research,  Development  and  Engineering  Center  (CRDEC).  The 
Computational  Aerodynamics  Branch,  Launch  and  Flight  Division,  U.S.  Army 
Ballistic  Research  Laboratory  (BRL)  has  been  tasked  to  determine  the  internal 
fluid  dynamics  within  a  partially,  liquid-filled  shell  subject  to  large  rapid¬ 
ly  applied  pressures  which  result  from  internal  burster  tube  rupture.  The 
main  area  of  concern  is  the  determination  of  the  time  history  of  the  internal 
pressure  on  the  shell  sidewalls  and  base. 


Previously,  an  analytical  technique  had  been  developed  for  predicting  the 
internal  fluid  dynamics  of  a  partially  liquid-filled,  nonspinning  cylinder 
subject  to  an  impulsively  applied  pressure  at  one  end.1  While  this  analytical 
technique  did  highlight  many  of  the  important  fluid  dynamic  processes  which 
occur  internally,  the  effect  of  important  parameters  such  as  the  internal 
shell  geometry  and  time  dependence  of  the  burster  tube  rupture  pressure  could 
not  be  considered. 


Due  to  the  difficulties  in  incorporating  these  parameters  into  an  analyt¬ 
ical  model,  a  numerical  model  has  been  developed  and  is  discussed  in  this  re¬ 
port.  This  numerical  model  solves  the  quasi  one-dimensional  Euler  equations 
and  a  equation  of  state  for  a  liquid  using  the  explicit  MacCormack  numerical 
technique.  Results  have  been  generated  for  an  internal  shell  geometry  cur¬ 
rently  being  tested  by  the  Terminal  Ballistics  Division  (TBD),  BRL.2  Compari¬ 
sons  are  made  between  the  numerically  predicted  pressure  levels  inside  the 
liquid  filled  shell  and  the  pressure  measurements  performed  by  TBD. 

II.  NUMERICAL  MODEL 

A  numerical  model  has  been  developed  to  simulate  the  internal  fluid  dyna¬ 
mics  of  a  nonspinning  partially  liquid-filled  shell  due  to  burster  tube  rup¬ 
ture.  In  this  section  the  numerical  model  is  presented.  First,  a  brief 
description  of  the  configuration  and  the  basic  fluid  dynamics  is  given;  the 
governing  equations  and  boundary  conditions  are  then  described;  and  finally 
the  numerical  procedure  is  outlined. 

1.  DESCRIPTION  OF  CONFIGURATION  AND  BASIC  FLUID  DYNAMICS 

Figure  1  depicts  the  internal  configuration  being  tested  by  the  Terminal 
Ballistics  Division  and  simulated  in  the  current  numerical  model.  The  inter¬ 
nal  configuration  is  an  ogive-cyl  inder-boattail  with  a  central  burster  tube 
which  extends  axially  from  the  nose.  The  axial  dimension  of  the  burster  tube 
is  approximately  half  the  internal  length  of  the  shell.  The  wall  of  the 
burster  tube  has  been  purposely  made  weaker  near  the  nose  to  control  the  loca¬ 
tion  of  burster  tube  rupture.  The  test  configuration  is  aligned  nose  down, 
leaving  some  percentage  of  ullage  (air  space  above  the  liquid)  in  the  boat- 
tailed  section  of  the  shell,  as  shown  in  Figure  2. 

The  fluid,  initially  at  rest,  is  set  in  motion  when  the  burster  tube  rup¬ 
tures  due  to  the  ignition  of  the  propellant  within  the  burster  tube.  The  rup¬ 
ture  of  the  burster  tube  subjects  the  surrounding  liquid  to  a  locally  high 
pressure  causing  a  pressure  wave  to  propagate  up  the  liquid,  as  shown  in 
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Figure  2.  This  pressure  wave  propagates  up  the  fluid  and  reflects  from  the 
liquid-air  interface.  After  the  reflection  of  the  pressure  wave  from  the 
liquid-air  interface,  the  interface  begins  moving  towards  the  shell  base.  A 
short  time  later  the  interface  impacts  on  the  base  and  another  pressure  wave 
is  generated  which  propagates  down  the  fluid  inside  the  shell. 

Several  features  of  the  fluid  motion  can  be  deduced  from  this  proposed 
description  of  the  fluid  motion;  (1)  the  compressible  nature  of  the  liquid  is 
important,  (2)  the  primary  motion  of  the  fluid  is  in  the  direction  along  the 
shell  axis  of  symmetry,  and  (3)  the  interfaces  which  bound  the  fluid  slug  move 
as  a  result  of  the  high  pressure  from  the  burster. 

2.  GOVERNING  EQUATIONS 

The  motion  of  the  fluid  within  the  shell  can  be  represented  by  the  time- 
dependent,  quasi  one-dimensional  Euler  equations,  shown  in  vector  form  below. 
This  vector  equation  represents  the  equations  of  mass  and  momentum  conserva¬ 
tion  along  the  axial  direction  within  the  shell;  where  p  is  the  local  fluid 
density,  u  is  the  axial  fluid  velocity,  and  P  is  the  local  fluid  pressure. 
The  internal  area  change  within  the  shell  is  accounted  for  through  the  source 

term,  ft,  where  A  is  the  local  cross  sectional  area  of  the  shell. 
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set  of  equations  is  accomplished  using  an  equation  of 
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where  E  is  the  bulk  modulus  of  the  fluid.  In  using  this  equation  of  state,  an 
isothermal  process  is  assumed.  A  fluid  bulk  modulus  of  320,000  psi  is  used 
throughout  the  computations.4 

The  equation  of  state  can  be  incorporated  within  the  equations  of  motion, 
yielding  the  following  vector  equation, 
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Because  this  set  of  equations  is  to  be  solved  on  a  grid  (x,t)  which  moves 
with  the  fluid  slug  as  a  function  of  time,  it  is  convenient  to  transform  these 
equations  into  a  computational  space  (£,  t)  which  remains  fixed  in  time  and 
space.  The  resulting  equations,  after  the  transformation  has  been  applied, 
are  written  below. 

(q)  +  (5t  q  +  Cx  F)  +  H  =  0 

where 

5ts  If  =  -  xt  5x  •  J  ■  1  /x5  (4) 


r  =  H  _  ,  „  =  3x  ..  =  3x 

s  -  Tx  ■ d  ’  xc  -  si  *  xt  ■  nr 


q=-q  ,  F  =  -  f  ,  H  =  I  ft  . 
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3.  BOUNDARY  CONDITIONS 

The  governing  equations  discussed  above  describe  the  fluid  motion  inside 
the  fluid  domain.  Using  the  quasi  one-dimensional  approach,  the  state  of  the 
fluid  must  also  be  specified  at  the  two  endpoints  of  the  domain. 

One  of  the  domain  endpoints  is  located  at  the  location  where  the  burster 
breaks.  At  this  point  the  fluid  is  exposed  to  high  pressure  from  the  burster. 
This  high  pressure  has  been  measured  experimentally  by  TBD  and  is  used  in  the 
numerical  model  to  specify  the  pressure  as  a  function  of  time  at  this  endpoint 
in  the  computational  domain.  Using  the  equation  of  state,  the  first  of  the 
two  dependent  variables,  the  density,  can  be  specified. 


p(t)  =  p0  exp  [P(t)/E]  (5) 

where  pQ  is  the  liquid  density  at  ambient  conditions.  The  second  dependent 
variable,  pu,  can  be  obtained  by  applying  the  momentum  equation. 

At  the  other  endpoint  of  the  domain,  the  liquid-air  interface,  one  of  two 
sets  of  boundary  conditions  is  applied  depending  of  the  state  of  the  inter¬ 
face.  The  first  set  of  boundary  conditions  is  applied  between  the  time  the 
burster  breaks  and  the  liquid-air  impacts  on  the  shell  base.  Here,  the 
pressure  acting  on  the  interface  is  nearly  ambient,  compared  with  the  pres¬ 
sures  produced  by  burster.  Using  the  equation  of  state,  the  first  of  the  two 
dependent  variables,  the  density,  can  be  specified, 


The  second  dependent  variable,  pu,  can  be  obtained  by  applying  the  momentum 
equation. 

The  second  set  of  boundary  conditions  is  applied  after  the  interface 
impacts  on  the  shell  base.  The  fluid  is  brought  to  rest  locally,  and  the 
second  dependent  variable  becomes  zero  at  this  endpoint  in  the  fluid  domain. 
The  boundary  condition  for  the  density  is  obtained  by  applying  the  continuity 
(mass  conservation)  equation. 

4.  NUMERICAL  PROCEDURE 

The  governing  equations  are  solved  using  the  second  order  accurate  explic¬ 
it  MacCormack  finite-difference  numerical  procedure.5  The  solution  is  obtain¬ 
ed  by  dividing  the  computational  domain  into  a  number  discrete  intervals  (or 
grid  points)  and  the  solution  advanced  in  time  using  the  following  predictor- 
corrector  approach. 
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where  superscript  "N"  denotes  the  time  level  of  the  solution  and  the  subscript 
"j"  denotes  the  spacial  location  of  the  grid  point. 

Boundary  conditions  are  updated  at  the  predictor  and  corrector  steps  using 
the  appropriate  one  sided  spacial  differencing.  After  the  solution  is 
updated,  the  updated  position  of  the  grid  is  determined  by  integrating  the 
velocity  of  the  grid  endpoints  in  time. 

The  solutions  presented  in  this  report  were  obtained  using  251  grid  points 
and  typically  required  1750  time  steps  to  span  the  time  interval  of  in¬ 
terest.  Each  solution  required  approximately  seven  minutes  of  CPU  time  on  a 
VAX  780  computer. 


III.  RESULTS 

The  results  of  the  application  of  the  numerical  model  are  presented  in  two 
sections.  The  first  section  is  concerned  with  validation  of  the  numerical 
model,  and  comparisons  are  made  between  the  numerical  model  and  experimental 
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results  obtained  by  TBD.  In  the  second  section,  results  of  a  parametric  study 
examining  the  effect  of  ullage  and  the  internal  shell  boattail  angle  are 
presented. 

1.  VALIDATION  OF  THE  NUMERICAL  MODEL 

A  total  of  18  cases  (various  ullages  and  burster  pressures  histories)  have 
been  run  using  the  numerical  model  and  comparisons  made  with  the  experimental 
pressure  data  obtained  by  TBD  to  provide  validation  of  the  numerical  model. 
Data  from  five  of  these  cases  are  presented  in  this  report.  TABLE  1  displays 
the  ullage  for  each  of  the  five  cases.  In  the  experiments  performed  by  TBD, 
the  driving  pressure  produced  by  the  burster  was  tailored  by  loading  the 
burster  with  varying  amounts  of  propellant.  The  amount  of  propellant  used  in 
each  test  is  also  shown  in  TABLE  1. 


TABLE  1.  Description  of  Cases. 


Case  ID 

Ul 1  age 

Burster 

FF3 

10% 

40gm 

FF4 

5% 

40gm 

FF5 

3% 

40gm 

FF7 

J  ,0 

50gi,i 

FF1 1 

3 /a 

30gm 

The  TBD  pressure  data  consists  of  output  from  five  pressure  gauges  located 
inside  the  pressure  vessel.  The  locations  of  four  of  these  gauges  is  shown  in 
Figure  3.  The  additional  gauge  (Gauge  1)  is  located  inside  the  burster  tube 
and  does  not  provide  information  of  interest  to  the  current  study.  The 
pressure  output  from  Gauge  2  is  used  to  specify  the  pressure  time-history  at 
one  endpoint  in  the  fluid  domain  in  the  numerical  model.  The  data  from  this 
gauge  is  numerically  filtered  before  each  computation  to  eliminate  high  fre¬ 
quency  oscillations  present  in  the  experimental  data.  (In  the  results  shown 
in  this  report  only  the  Gauge  2  experimental  data  has  been  smoothed.  Experi¬ 
mental  data  shown  from  the  other  gauges  is  the  raw  experimental  data.) 

Case  FF5  is  examined  first  in  some  detail  to  reconcile  the  physics  with 
the  experimental  and  numerical  output.  Results  from  the  other  cases  follow  a 
similar  pattern.  Figure  4a  displays  the  experimental  output  of  Gauge  2  and 
the  smoothed  data  used  as  input  to  the  numerical  model.  Figures  4b-4d  show 
the  experimentally  measured  and  numerically  predicted  pressure  at  Gauges  3,  4 
and  5  respecti vely ,  and  reference  to  these  figures  is  made  repeatedly  in  the 
following  several  paragraphs.  In  order  to  aid  in  understanding  the  results  in 
these  figures,  critical  events  have  been  marked  on  the  graphs  by  alphabetic 
characters,  and  are  referenced  to  in  the  text. 

The  event  begins  with  the  breakage  of  the  burster  tube.  Gauge  2,  which  is 
near  the  location  of  burster  rupture,  experiences  a  rapid  increase  in  pressure 
almost  immediately.  Two  tenths  of  a  millisecond  later  the  primary  pressure 
wave  has  propagated  up  the  pressure  vessel  to  Gauge  3(A),  and  at  0.35  milli¬ 
seconds  the  pressure  wave  passes  Gauge  4(B).  The  pressure  wave  continues  to 
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show  that  between  0.7  and  0.85  milliseconds  the  lower  boundary  of  the  fluid 
domain  has  been  pushed  past  Gauge  3  by  the  burster  gases.  During  this  time 
period,  high  frequency  oscillations  in  the  experimental  pressure  trace  are 
evident,  which  may  suggest  that  experimentally  the  gauge  is  being  exposed  to 
the  combustion  gases  from  the  burster  tube.  The  computed  and  experimental 
pressure  history  at  Gauge  4  are  in  fairly  good  agreement.  The  numerically 
predicted  pressures  at  Gauge  5  are  slightly  greater  than  the  experimental 
values  and  again  the  pressure  relief  on  the  shell  base  occurs  later  in  the 
computation  than  in  the  experiment. 

The  pressure  traces  for  Case  FF11  are  shown  in  Figures  8a,b,c.  The  com¬ 
puted  pressures  at  Gauges  3  and  4  are  again  in  good  agreement  with  the  exper¬ 
imental  values.  The  computed  pressures  at  Gauge  5  are  significantly  higher 
than  the  experimentally  measured  values,  though  the  experimental  results 
indicate  negative  pressures  after  1.0  millisecond  indicating  a  possible  gauge 
zero-drift  problem. 

In  general,  the  numerical  results  show  excellent  agreement  with  the  exper¬ 
imental  results  at  Gauge  3.  This  is  not  surprising  due  to  the  close  proximity 
of  Gauge  3  to  Gauge  2,  which  is  used  as  the  driving  pressure  for  the  numerical 
model.  The  numerical  results  indicate  that  in  a  number  of  cases  the  lower 
liquid  interface  is  pushed  past  gauge  location  3,  particularly  as  the  ullage 
is  increased.  The  numerical  predictions  of  the  pressure  history  at  Gauge  4 
are  in  good  agreement  with  the  experimental  measurements.  The  numerical 
results  at  Gauge  5  display  qualitatively  the  same  features  as  does  the  experi¬ 
mental  results.  The  arrival  time  of  the  fluid  interface  is  well  predicted  and 
in  many  cases  the  pressures  levels  due  to  the  impacting  of  the  fluid  interface 
on  the  base  agree  well  with  the  experimental  results.  The  eventual  relief  in 
pressure  on  the  shell  base  occurs  somewhat  later  in  the  computation  than  in 
the  experiment,  although  this  does  not  appear  to  be  a  critical  parameter. 

2.  PARAMETRIC  STUDY 

One  of  the  primary  advantages  of  developing  a  numerical  model,  such  as  the 
one  discussed  here,  is  that  parametric  studies  can  be  quickly  and  easily 
performed.  The  code  developed  in  this  study  has  been  used  to  elicit  the  im¬ 
portant  parameters  of  the  forced  liquid  problem.  Two  parameters  examined  here 
are  the  effects  of  ullage  and  boattail  angle. 

a.  Ullage  A  parametric  study  examining  the  effect  of  ullage  on  the 
internal  pressure  history  has  been  performed.  Results  have  been  generated  for 
a  single  input  pressure  history  (corresponding  to  Case  FF5)  and  for  several 
ullages;  3,  5,  10  and  15  percent. 

Figures  9a  and  9b  show  the  pressure  history  at  Gauge  3  for  5  and  10 
percent  ullages.  For  the  first  0.6  milliseconds,  the  pressure  history  shows 
little  variation  between  the  two  ullages.  This  was  also  observed  in  the 
computed  results  for  3  and  15  percent  ullages.  After  0.6  milliseconds,  the 
numerical  results  show  that  the  lower  pressurized  fluid  surface  is  pushed  past 
Gauge  3  for  a  duration  of  time  that  increases  as  the  ullage  increases.  This 
behavior  was  observed  in  the  computed  results  across  the  range  of  ullages 
exami ned. 
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The  pressure  histories  at  Gauge  5  for  each  of  the  ullages  are  shown  in 
Figures  10a,b,c,d.  The  figures  show  that  the  base  is  generally  exposed  to 
higher  pressures  as  the  ullage  is  increased.  However,  the  base  experiences 
the  pressure  loading  over  a  smaller  duration  as  the  ullage  is  increased. 

b.  Internal  boattail  angle  A  study  was  performed  to  demonstrate  the 
effect  of  the  internal  boattai 1  angle  of  the  shell  on  the  internal  pressure 
history.  Computations  were  performed  for  two  boattail  angles  (7  and  10 
degrees)  using  the  FF5  input  pressure,  and  a  ullage  of  10%. 

The  pressure  history  at  Gauge  3  show  very  little  variation  due  to  the 
change  in  the  internal  boattail  angle.  The  pressure  history  at  Gauge  4,  shown 
in  Figure  11a,  also  show  only  small  variations.  The  largest  effect  is  noted 
at  gauge  location  5  (on  the  shell  base)  where  larger  peak  pressures  are  pre¬ 
dicted  for  the  larger  boattail  angle,  as  shown  in  Figure  lib.  The  duration 
and  average  magnitude  of  the  pressure  loading  are  approximately  the  same. 


IV.  CONCLUSIONS 

A  quasi  one-dimensional  numerical  model  has  been  developed  to  compute  the 
internal  fluid  dynamics  of  a  nonspinning  partially  liquid-filled  shell  which 
results  from  burster  tube  rupture.  Computational  results  have  been  obtained 
for  a  shell  configuration  currently  being  tested  by  TBD.  The  results  show 
very  good  agreement  when  compared  with  pressure  measurements  obtained  by  TBD, 
validating  the  numerical  approach.  The  numerical  model  has  provided  further 
insight  into  the  physical  processes  which  occur  within  the  shell. 

A  parametric  study  has  been  performed  using  the  numerical  model  to  examine 
the  effect  of  ullage  and  internal  boattail  angle  of  the  shell  on  the  internal 
surface  pressures.  The  numerical  results  show  that  the  fluid  interface 
located  in  the  vicinity  of  the  burster  tube  rupture  is  pushed  further  from  the 
nose  of  the  shell  by  the  high  pressure  gases  as  the  ullage  is  increased.  The 
results  also  show  that  the  shell  base  is  generally  exposed  to  higher  pressures 
as  the  ullage  is  increased,  however,  the  base  experiences  these  higher  pres¬ 
sures  over  a  somewhat  shorter  duration.  The  numerical  results  show  that  as 
the  boattail  angle  of  the  shell  is  increased  the  peak  base  pressures  are 
increased,  though  the  average  pressure  and  the  duration  of  the  pressure  load 
on  the  shell  base  remain  constant. 
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LIST  OF  SYMBOLS 


A  =  local  cross  sectional  area 

E  =  fluid  bulk  modulus 

F  =  flux  vector 

F  =  flux  vector  normalized  by  Jacobian 

H  =  source  term  due  to  area  change 

H  =  source  term  normalized  by  Jacobian 

J  =  Jacobian  of  the  transformation 

P  =  pressure 

q  =  vector  of  dependent  variables 

q  =  vector  of  dependent  variables  scaled  by  Jacobian 

t  =  time 

u  =  axial  velocity 

x  =  axial  coordinate 

p  =  density 

pQ  =  density  at  ambient  conditions 

Subscri pts 

t  =  derivative  with  respect  to  time 

x  =  derivative  with  respect  to  axial  coordinate 

£  =  derivative  with  respect  to  spacial  computational  coordinate 

Superscripts 

N  =  denotes  solution  at  "N"th  time  level 

N+l  =  denotes  solution  at  "N+T'tb  time  level 

N+l  =  denotes  solution  at  end  of  "N+l"th  predictor  step 
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